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Abstract 

We present a stochastic model which describes fronts of cells in- 
vading a wound. In the model cells can move, proliferate, and ex- 
perience cell-cell adhesion. We find several qualitatively different 
regimes of front motion and analyze the transitions between them. 
Above a critical value of adhesion and for small proliferation large 
isolated clusters are formed ahead of the front. This is mapped onto 
the well-known ferromagnetic phase transition in the Ising model. 
For large adhesion, and larger proliferation the clusters become con- 
nected (at some fixed time). For adhesion below the critical value the 
results are similar to our previous work which neglected adhesion. 
The results are compared with experiments, and possible directions 
of future work are proposed. 



1 Introduction 

When a wound heals, surrounding cells fill the wounded area by enhanced 
motility (i.e. diffusion) and enhanced proliferation (see Ref. for a 

recent review) . Most theoretical treatments of this process [SJ 1141 ITS] 
employ a reaction-diffusion equation for the cell density, equivalent to the 
Fisher-Kolmogorov (FK) equation [7] . Another approach was taken in 
m^], where a very simple discrete model was formulated for the similar 
problem of flame-front propagation. Recently, this model was applied to 
wound healing 2 . The model takes into account proliferation and diffusion, 
and for small proliferation it reduces to the FK equation. Biologically 
reasonable proliferation rates are small compared to rates of diffusion 
so the front velocity is in a good agreement with experimental findings both 
for continuum and discrete models. However, the theoretically predicted 



1 



width of the front is much larger than the one measured experimentally 
(see, for example, Ref. (17)1. 

Walker et al. |171 ITS] have proposed that the answer to the paradox 
lies in the inclusion of cell-cell adhesion. They investigated an agent based 
model and observed qualitatively different regimes of cell organization for 
low and high values of adhesion. In order to investigate this idea further, 
we consider a simple discrete model which describes the phenomenon of 
wound healing, focusing on the key processes: cell-cell adhesion, diffusion, 
proliferation. Simulations of our model show two qualitatively different 
regimes depending on the adhesion strength, q, similarly to the results of 
Walker et al. ^1 ■ We found that the transition between these regimes 
is very sharp, and related to the phase transition in the Ising model (see 
below). Another regime of cell organization was found, depending on the 
proliferation rate, a. Here we report preliminary results of our study and 
analyze the transitions between the regimes. Finally we discuss the biolog- 
ical applications of our work and compare the results with experiments. 

2 Formulation of the model 

We will formulate our model in a way which is reminiscent of the standard 
'scratch assay' experiment in wound healing studies [T5|. Consider a square 
two-dimensional lattice in a channel geometry. Each lattice site can be 
empty or once occupied by a cell. We assume the lattice distance to be 
equal to cell diameter (of the order of 10/im), taking into account hard-core 
exclusion. Thus, a fully occupied region of the lattice represents a confluent 
monolayer, that is, unwounded or healed tissue. 

Initially, we put cells into the left part of the channel. We take x 
to measure distance along the channel. In the initial state all sites with 
x < 40 are occupied and the rest empty. Thus x — 40 is the edge of the 
wound. For t > cells diffuse and proliferate along the channel. The 
dynamical rules which define the model are as follows: A cell is picked at 
random, and one of the four neighboring sites is also picked at random. 
If this site is empty, the cell can proliferate to this site (so a new cell is 
born there), or migrate there. The probability for proliferation is a. The 
probability for migration decreases with the number of nearest neighbors 
so that p m i gr — (1 — a)(l — q) n , where < q < 1 is the adhesion parameter, 
and 1 < n < 4 is the number of nearest neighbors. The case q = means 
no adhesion and brings us back to the simple model For nonzero 

q, it is much harder to a cell to diffuse if it has many neighbors. After each 
step time is advanced by 1/iV, where N is the current number of cells. 

As healing proceeds there is a zone in front of the healed tissue which 
is partly filled with cells. We call this the invasive region. We have done 
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extensive simulations of this model. We report the results in the following 
sections. 



3 Structure of the invasive region 

Our simulations show that the dynamics and structure of the invasive zone 
is qualitatively different depending on the two parameters of the model: a 
and q. Figure^shows the (a, q) phase plane and points out three different 
regions of parameters. 
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Figure 1: Phase plane (a, q). A qualitatively different behavior is observed 
in different regions (a,b,c) in the phase plane, see Fig. [2 Two arrows denote 
transitions between (a) and (b), and between (b) and (c), see text. 

The different types of behavior are shown in Figure[5]by means of three 
snapshots of the system which correspond to points (a), (b), and (c) in 
Figure ^ Figure shows the system for small proliferation and small 
(subcritical, see below) adhesion, region (a) in Fig. ^ Here, a front of 
cells propagates along the channel, both front velocity and front width are 
well-defined (and can be calculated as in 0). Figure |5Jd shows a snapshot 
for the same proliferation rate and large (supercritical, see below) adhesion 
strength, region (b) in Fig. ^ fn contrast to the propagating fronts shown 
in Fig. one can see a number of isolated clusters that are formed in 
the invasive region. However, as we increase proliferation (for the fixed 
adhesion parameter), another interesting transition occurs. For moderate 
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proliferation, large clusters in the invasive zone become connected to each 
other and to the initial dense front, as can be seen in Figure^, region (c) 
in Fig. [T] 
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Figure 2: Snapshots of the system in three qualitatively different regions in 
the phase plane of parameters, see Fig. 2] Each white dot corresponds to an 
occupied site and each black dot to an empty site; healing proceeds from left 
to right. The parameters are q = 0.7, a = 3x 10 -4 (a), q = 0.9, a = 3x 10~ 4 
(b), and q = 0.9, a = 7 x 10~ 4 (c). 



We can qualitatively analyze the two transitions which are shown by 
arrows in Fig. ^ First we focus on the transition from (a) to (b) which 
occurs at a fixed (and sufficiently small, see below) proliferation, when the 
adhesion parameter crosses critical value. 

We point out that our model without proliferation can be mapped into 
the Ising model in statistical physics. In this mapping, an empty site corre- 
sponds to spin "down", an occupied site corresponds to spin "up", so that 
there is a simple relation between the average density, c, and the average 
magnetization, m, in the Ising model: c = (m+l)/2. The adhesion param- 
eter q is related to the ratio of magnetic coupling J and the temperature T 
by q = 1 — exp(— J/ksT). The mapping is possible because our dynamical 
rules satisfy detailed balance. Therefore, the statics of our model is the 
same as in the Ising model. By statics, we mean a phase diagram (m, T) 
(or (c, q) in our case) which has stable and unstable regions. In the stable 
region, a homogenous state (with uniformly distributed cells) remains ho- 
mogenous; in contrast, in the unstable region phase separation occurs and 
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Figure 3: The critical adhesion parameter as a function of density as given 
by Eq. (I). An inset shows the derivative dq/dc versus c. 

The two-dimensional Ising model was solved by Onsager |1 1|. and the 
curve m(T), which separates the stable and unstable regions, is known. In 
terms of average density c and the critical adhesion parameter q c , we have: 



The unstable region corresponds to q > g c , so for supercritical adhesion 
large clusters are formed. Of course, our system is not homogenous. How- 
ever the density dependence of q c is rather slow in a wide range of interme- 
diate densities, see Figure|3] Therefore, we can roughly estimate the critical 
adhesion as q c ~ q c (c — 0.5) = 0.8284. This is the horizontal dashed line in 
Figure^] This observation explains the different structures observed in the 
invasive zones of Figs. [2k andUJa (the transition from region (a) to region 
(b), Figured). 

The propagation of fronts is a very important topic in statistical me- 
chanics. It turns out that there are propagating fronts in region (a), similar 
to those of the FK equation |7J. We averaged over a series of simula- 
tions (and over the channel width) to obtain smooth density profiles in 
different parameter regions. Figure 0] shows an example of such a front 
and a velocity of the fronts as a function of adhesion parameter for dif- 
ferent values of proliferation. Note a rather slow velocity dependence on 
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the adhesion parameter q. Clearly, as q goes to zero, the front velocity v 
approaches its theoretical value from the FK equation, which is given in 
our notations by a 1 / 2 . Another important issue is the scaling properties of 
the front. For example, it would be interesting to analyze whether there 
is a KPZ roughening as in the case of zero adhesion Unfortunately, 
it is extremely difficult to investigate fronts roughening in this problem, 
because the proliferation is very small, so an intrinsic front width is quite 
large (see Figure 0}. 
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Figure 4: Front velocity as a function of adhesion parameter for different 
values of proliferation. The dotted curve corresponds to a = 7 x 10~ 4 , the 
dashed curve corresponds to a = 3 x 10~ 4 . An inset shows a typical density 
profile in region (a). Here q = 0.8, a = 3 x 10~ 4 , and the time is t = 10 5 . 



Now we turn to the transition from region (b) to region (c), see Fig.^ 
We increase the proliferation parameter a, keeping the adhesion parameter 
q > q c constant. In order to compare the results for different values of 
proliferation, we fix the total time of the simulations, t = 3 x 10 4 , and 
measure the maximum distance L (in x direction) one can move through 
occupied sites. In other words, consider a path that passes only through 
occupied sites. L is just the x component of the longest path along the 
channel. Figure [5] shows the dependence of L on a. For small proliferation 
rate, L is small (as the clusters in the invasive zone are not connected to 
each other and to initial dense region x < 40) and it slowly increases with 
a up to a transition point. But then a small increase in a is followed by a 
rapid increase in L as large clusters become connected, see Figs. |2t- Note 
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Figure 5: Length of the maximal path along the channel as a function of 
proliferation, see text. Note that the proliferation threshold depends on 
time: a larger proliferation is needed to get this percolation-like transition 
at a smaller time. The adhesion parameter is q = 0.9, the time is t = 3x 10 4 . 



that the proliferation threshold depends on time: larger proliferations are 
needed to get this percolation-like transition at smaller times. Therefore, 
the diagonal dashed line in Fig. ^ which determines the transition, is not 
fixed and moves to the right (to larger proliferations) for smaller times. 
This indicates that region (b) is a transient region. For any (biologically 
reasonable) proliferation and for supercritical adhesion the system is in 
region (b) at early times. However, at very late times there are propagating 
fronts with some defined width and velocity. The fronts develop well inside 
region (c), much later than the transition from (b) to (c) occurs. So, 
the transition from (b) to (c) is sharp but occurs in a transient regime. 
The connection to percolation problem would be an interesting direction 
of future research. 

To emphasize the effect of cell-cell adhesion, we performed simulations 
in the geometry of a scratch assay experiment |17| . (Up until now we have 
considered only one side of the wound.) For the type of cells used by Walker 
et al. ^| adhesion is controlled by the concentration of Ca ++ in the 
system. Lower concentrations of calcium suppress adhesion, while higher 
concentrations promote it. 

In the experiments |17| an initial width of a scratch made in a mono- 
layer of epithelial cells was about 600/im width, which roughly corresponds 
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Figure 6: Simulations of wound healing. The upper panel corresponds to 
small adhesion q = 0.7, the lower panel to larger adhesion q = 0.99. White 
dots represent occupied sites and black dots empty sites. The dashed lines 
indicate initial wound edges. The proliferation parameter is a = 0.01. 

to 60 cell diameters. Since the characteristic proliferation time (of the or- 
der of 1 day) is much larger than the characteristic hopping time (of the 
order of minutes), the basic time unit in our model approximately equals 
the hopping time. Therefore we ran our model out to a time of 10 3 , which 
corresponds to about 24 hours in the experiment ^7]. Figure El shows the 
results of simulations both for small adhesion strength (upper panel) and 
large adhesion strength (lower panel) for the same value of proliferation 
rate. For small adhesion the front is wide and wound closure is almost 
complete, whereas for large adhesion there is a very sharp and slowly mov- 
ing front, in a good agreement with the experiments (Figs. 7 and 8 in 
Ref. |HJ) both for high and low adhesion. There is another point that we 
would like to clarify. Figure shows clusters of cells separated by cell- free 
patches, while Figure (lower panel) does not. This is due to the fact that 
one should consider a much larger system (by factor of 6) and wait for a 
much longer time (by factor of 30) to get the structure shown in Figure 
13;. 

4 Summary and discussion 

In this work we presented a stochastic discrete model of wound healing. 
In contrast to previous modeling attempts, both discrete and continuum, 



8 



we took into account the phenomenon of cell-cell adhesion. Large adhesion 
gives rise to sharp wound healing fronts even for very low proliferation, 
see Fig. [S] In general, depending on the two parameters of the model, the 
adhesion q and the proliferation a, a completely different behavior in the 
invasive zone can be observed. 

Certainly, our model oversimplifies a complex biological process of wound 
repair. There are many complex biochemical processes such as calcium 
waves, apoptosis, etc. which we do not attempt to treat. One feature 
which we can treat is the fact that the model does not take into account a 
real biological cell cycle which may influence the dynamics. A cell cycle can 
be incorporated into the model in the following way which is based on the 
work of Walker et al. ^| E| • Each cell is treated as an agent with its own 
internal time. A cell moves and adheres according to the rules described 
above. Additionally, its internal clock advances when it is sampled. When 
a cell reaches the end of the cell cycle length, mitosis occurs and the cell 
ceases to move for a small fraction of the cycle time. It then proliferates 
in a random direction with a given probability. In the limit of long times, 
the parameter a used in the simple model corresponds to the product of 
the division probability and the cycle frequency. With randomized initial 
internal times, this model behaves very similarly to the purely probabilistic 
approach which we have discussed above. 

Now we would like to emphasize our main new results, comparing to 
Refs. |17l I18| . One important new finding, resulting from our theoretical 
analysis, is the existence of a sharp transition between qualitatively differ- 
ent types of behavior in cell cultures when an adhesion parameter passes 
a certain value. Though qualitatively different structures for low and high 
adhesion were observed earlier (Refs. the fact that a very small 

change in adhesion parameter may give a qualitatively different behavior 
was unknown. This result shows a deep and surprising analogy between 
ferromagnetic phase transition in statistical physics and transition to clus- 
tering in ensemble of cells. This analogy had not also been found prior to 
our work. 

Another new result is the transition from region (b) to region (c), which 
occurs when increasing proliferation for supercritical adhesion. This tran- 
sition and the detailed clusters structure for supercritical adhesion had not 
been investigated in Refs. |17l IIS] . Finally, there is another new result, 
showing a typical front profile of cells density and the front velocity as a 
function of adhesion parameter for different values of proliferation. That 
possibly gives a way to theoretical treatment of the problem of cell invasion 
into a wound in terms of a modified Cahn-Hilliard equation. 

Indeed, a promising direction of future work is continuum modeling of 
front propagation, both for subcritical and supercritical adhesion. A proper 
candidate here is Cahn-Hilliard equation, which describes the dynamics of 
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phase separation below the critical temperature ^JEIE]- hi our system 
this corresponds to supercritical adhesion and zero proliferation. Indeed, 
one can show jSJ^l that a variant of Cahn-Hilliard equation can be derived 
from discrete lattice gas model with nearest neighbors interaction. This 
equation applies to models with conserved order parameter (in our case, 
the model without proliferation). One can also try to take into account 
proliferation similarly to the FK equation, adding term ac(l — c) to the 
Cahn-Hilliard equation. This work is in progress. 
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